brms today! Install this package
now, if you haven’t already.install.packages(c("brms","tidybayes"))
library(here)
library(tidyverse)
library(cowplot)
library(brms)
library(tidybayes)
library(patchwork)
# here's a bunch of stuff I include at the beginning of each qmd file.
# it's mostly setting the color palette for figures, but it also controls the default code chunk option (showing you, i.e., echo) and setting the figure quality to ensure that it looks nice
knitr::opts_chunk$set(fig.retina=3, echo=TRUE)
theme_set(theme_cowplot())
default_palettes <- list(
c("#5e8485" , "#0f393a") ,
c("#1c5253" , "#5e8485" , "#0f393a") ,
# palette with 5 colours
c( "#1c5253" , "#e07a5f", "#f2cc8f" , "#81b29a" , "#3d405b" ) ,
# same palette interpolated to 8 colours
c( "#1c5253" , "#e07a5f", "#f2cc8f" , "#81b29a" , "#3d405b" , "#a7a844" , "#69306d" )
)
options(ggplot2.discrete.fill = default_palettes,
ggplot2.discrete.colour = default_palettes)
State a clear question.
Sketch your causal assumptions.
Use the sketch to define a generative model.
Use the generative model to build an estimator.
Profit.
Here’s an example:
\[\begin{align*} y_i &\sim \text{Normal}(\mu_i,\sigma) \\ \mu_i &= \beta x_i \\ \beta &\sim \text{Normal}(0,10) \\ \sigma &\sim \text{Exponential}(1) \\ x_i &\sim \text{Normal}(0,1) \\ \end{align*}\]
Here’s the model for last week’s globe-tossing experiment:
\[\begin{align*} W &\sim \text{Binomial}(N,p) \\ p &\sim \text{Uniform}(0,1) \\ \end{align*}\]
The whole model can be read as:
The count \(W\) is distributed binomially with sample size \(N\) and probability \(p\). The prior for \(p\) is assumed to be uniform between 0 and 1.
Here’s the model for last week’s globe-tossing experiment:
\[\begin{align*} W &\sim \text{Binomial}(N,p) \\ p &\sim \text{Uniform}(0,1) \\ \end{align*}\]
brmsLast week, we used grid approximation to estimate the posterior distribution. In the video you watched for today, McElreath moves on to something called QUADRATIC APPROXIMATION. It’s good to understand what that’s doing, but you and I are moving right along to MCMC. We won’t go into the details of what’s happening for a few weeks, but let’s start with the code to estimate the model.
Let’s say we tossed the globe 9 times and observed 6 waters:
m1 <-
brm(data = list(w = 6), # <1>
family = binomial(link = "identity"), # <2>
w | trials(9) ~ 0 + Intercept, # <3>
prior(uniform(0, 1), class = b), # <4>
iter = 5000, warmup = 1000, seed = 3, chains=1, # <5>
file = here("files/models/m21.1")) # <6>
Grid approximation gave us the calculated probability of each possible value of our parameter, \(p\). But our method of conducting bayes will no longer give us such a neat solution. Here’s how you get the posterior distribution for \(p\):
samples_from_post = as_draws_df(m1)
samples_from_post
## # A draws_df: 4000 iterations, 1 chains, and 3 variables
## b_Intercept lprior lp__
## 1 0.46 0 -2.1
## 2 0.49 0 -1.9
## 3 0.56 0 -1.5
## 4 0.66 0 -1.3
## 5 0.71 0 -1.3
## 6 0.62 0 -1.3
## 7 0.68 0 -1.3
## 8 0.55 0 -1.5
## 9 0.63 0 -1.3
## 10 0.67 0 -1.3
## # ... with 3990 more draws
## # ... hidden reserved variables {'.chain', '.iteration', '.draw'}
samples_from_post %>%
ggplot(aes(x=b_Intercept)) +
geom_density(fill = "grey", color = "white") +
labs(x="Proportion water")
ppd = posterior_predict(m1)
dim(ppd)
## [1] 4000 1
ppd
## [,1]
## [1,] 6
## [2,] 4
## [3,] 6
## [4,] 5
## [5,] 6
## [6,] 6
## [7,] 6
## [8,] 4
## [9,] 7
## [10,] 7
## [11,] 5
## [12,] 7
## [13,] 8
## [14,] 7
## [15,] 8
## [16,] 9
## [17,] 7
## [18,] 9
## [19,] 8
## [20,] 9
## [21,] 7
## [22,] 9
## [23,] 8
## [24,] 7
## [25,] 6
## [26,] 6
## [27,] 3
## [28,] 4
## [29,] 5
## [30,] 3
## [31,] 5
## [32,] 8
## [33,] 4
## [34,] 7
## [35,] 3
## [36,] 2
## [37,] 5
## [38,] 5
## [39,] 5
## [40,] 5
## [41,] 6
## [42,] 8
## [43,] 7
## [44,] 6
## [45,] 6
## [46,] 6
## [47,] 7
## [48,] 6
## [49,] 7
## [50,] 9
## [51,] 7
## [52,] 3
## [53,] 6
## [54,] 7
## [55,] 4
## [56,] 7
## [57,] 5
## [58,] 2
## [59,] 6
## [60,] 6
## [61,] 7
## [62,] 6
## [63,] 4
## [64,] 6
## [65,] 6
## [66,] 6
## [67,] 8
## [68,] 7
## [69,] 6
## [70,] 4
## [71,] 6
## [72,] 5
## [73,] 8
## [74,] 5
## [75,] 5
## [76,] 3
## [77,] 2
## [78,] 6
## [79,] 9
## [80,] 5
## [81,] 6
## [82,] 7
## [83,] 5
## [84,] 5
## [85,] 6
## [86,] 7
## [87,] 5
## [88,] 5
## [89,] 6
## [90,] 6
## [91,] 6
## [92,] 3
## [93,] 5
## [94,] 5
## [95,] 5
## [96,] 6
## [97,] 7
## [98,] 7
## [99,] 7
## [100,] 8
## [101,] 5
## [102,] 7
## [103,] 8
## [104,] 6
## [105,] 6
## [106,] 8
## [107,] 5
## [108,] 6
## [109,] 6
## [110,] 3
## [111,] 5
## [112,] 4
## [113,] 5
## [114,] 8
## [115,] 9
## [116,] 3
## [117,] 6
## [118,] 5
## [119,] 7
## [120,] 7
## [121,] 8
## [122,] 3
## [123,] 5
## [124,] 7
## [125,] 7
## [126,] 7
## [127,] 5
## [128,] 6
## [129,] 5
## [130,] 6
## [131,] 9
## [132,] 8
## [133,] 4
## [134,] 5
## [135,] 4
## [136,] 8
## [137,] 4
## [138,] 7
## [139,] 7
## [140,] 5
## [141,] 9
## [142,] 5
## [143,] 6
## [144,] 5
## [145,] 6
## [146,] 8
## [147,] 6
## [148,] 6
## [149,] 5
## [150,] 7
## [151,] 6
## [152,] 9
## [153,] 6
## [154,] 7
## [155,] 8
## [156,] 4
## [157,] 6
## [158,] 7
## [159,] 5
## [160,] 3
## [161,] 6
## [162,] 3
## [163,] 5
## [164,] 8
## [165,] 6
## [166,] 6
## [167,] 3
## [168,] 6
## [169,] 7
## [170,] 9
## [171,] 5
## [172,] 5
## [173,] 5
## [174,] 8
## [175,] 7
## [176,] 3
## [177,] 9
## [178,] 9
## [179,] 6
## [180,] 7
## [181,] 5
## [182,] 7
## [183,] 9
## [184,] 7
## [185,] 6
## [186,] 5
## [187,] 8
## [188,] 6
## [189,] 6
## [190,] 5
## [191,] 8
## [192,] 8
## [193,] 6
## [194,] 7
## [195,] 4
## [196,] 5
## [197,] 6
## [198,] 7
## [199,] 8
## [200,] 8
## [201,] 9
## [202,] 4
## [203,] 8
## [204,] 7
## [205,] 3
## [206,] 6
## [207,] 4
## [208,] 6
## [209,] 7
## [210,] 7
## [211,] 7
## [212,] 4
## [213,] 9
## [214,] 5
## [215,] 8
## [216,] 8
## [217,] 9
## [218,] 6
## [219,] 7
## [220,] 5
## [221,] 6
## [222,] 4
## [223,] 3
## [224,] 6
## [225,] 4
## [226,] 3
## [227,] 4
## [228,] 9
## [229,] 6
## [230,] 6
## [231,] 8
## [232,] 7
## [233,] 5
## [234,] 5
## [235,] 7
## [236,] 6
## [237,] 9
## [238,] 7
## [239,] 5
## [240,] 3
## [241,] 3
## [242,] 6
## [243,] 7
## [244,] 8
## [245,] 6
## [246,] 6
## [247,] 1
## [248,] 8
## [249,] 4
## [250,] 5
## [251,] 6
## [252,] 8
## [253,] 7
## [254,] 9
## [255,] 8
## [256,] 7
## [257,] 4
## [258,] 7
## [259,] 7
## [260,] 4
## [261,] 7
## [262,] 8
## [263,] 5
## [264,] 8
## [265,] 4
## [266,] 8
## [267,] 5
## [268,] 7
## [269,] 6
## [270,] 8
## [271,] 6
## [272,] 6
## [273,] 9
## [274,] 5
## [275,] 5
## [276,] 4
## [277,] 9
## [278,] 8
## [279,] 7
## [280,] 7
## [281,] 6
## [282,] 3
## [283,] 5
## [284,] 1
## [285,] 2
## [286,] 1
## [287,] 1
## [288,] 9
## [289,] 6
## [290,] 8
## [291,] 7
## [292,] 9
## [293,] 9
## [294,] 7
## [295,] 4
## [296,] 6
## [297,] 8
## [298,] 6
## [299,] 1
## [300,] 2
## [301,] 5
## [302,] 7
## [303,] 6
## [304,] 7
## [305,] 7
## [306,] 6
## [307,] 8
## [308,] 7
## [309,] 6
## [310,] 5
## [311,] 4
## [312,] 7
## [313,] 7
## [314,] 4
## [315,] 6
## [316,] 8
## [317,] 5
## [318,] 7
## [319,] 5
## [320,] 4
## [321,] 4
## [322,] 4
## [323,] 4
## [324,] 3
## [325,] 5
## [326,] 7
## [327,] 5
## [328,] 5
## [329,] 4
## [330,] 8
## [331,] 4
## [332,] 7
## [333,] 8
## [334,] 5
## [335,] 8
## [336,] 8
## [337,] 8
## [338,] 5
## [339,] 2
## [340,] 3
## [341,] 5
## [342,] 6
## [343,] 7
## [344,] 7
## [345,] 7
## [346,] 4
## [347,] 8
## [348,] 9
## [349,] 8
## [350,] 7
## [351,] 8
## [352,] 6
## [353,] 4
## [354,] 4
## [355,] 8
## [356,] 7
## [357,] 3
## [358,] 1
## [359,] 3
## [360,] 6
## [361,] 4
## [362,] 4
## [363,] 7
## [364,] 4
## [365,] 6
## [366,] 4
## [367,] 7
## [368,] 4
## [369,] 9
## [370,] 8
## [371,] 6
## [372,] 3
## [373,] 3
## [374,] 3
## [375,] 2
## [376,] 3
## [377,] 4
## [378,] 8
## [379,] 2
## [380,] 6
## [381,] 9
## [382,] 6
## [383,] 9
## [384,] 7
## [385,] 7
## [386,] 7
## [387,] 6
## [388,] 5
## [389,] 3
## [390,] 7
## [391,] 7
## [392,] 7
## [393,] 7
## [394,] 4
## [395,] 8
## [396,] 9
## [397,] 5
## [398,] 5
## [399,] 8
## [400,] 7
## [401,] 6
## [402,] 6
## [403,] 7
## [404,] 8
## [405,] 4
## [406,] 6
## [407,] 6
## [408,] 7
## [409,] 9
## [410,] 8
## [411,] 8
## [412,] 7
## [413,] 6
## [414,] 4
## [415,] 7
## [416,] 6
## [417,] 4
## [418,] 4
## [419,] 7
## [420,] 4
## [421,] 5
## [422,] 7
## [423,] 2
## [424,] 6
## [425,] 8
## [426,] 6
## [427,] 3
## [428,] 8
## [429,] 6
## [430,] 5
## [431,] 6
## [432,] 3
## [433,] 4
## [434,] 7
## [435,] 7
## [436,] 8
## [437,] 7
## [438,] 3
## [439,] 3
## [440,] 8
## [441,] 6
## [442,] 8
## [443,] 7
## [444,] 7
## [445,] 7
## [446,] 4
## [447,] 5
## [448,] 5
## [449,] 4
## [450,] 7
## [451,] 3
## [452,] 5
## [453,] 6
## [454,] 2
## [455,] 6
## [456,] 3
## [457,] 8
## [458,] 7
## [459,] 9
## [460,] 6
## [461,] 8
## [462,] 7
## [463,] 8
## [464,] 5
## [465,] 8
## [466,] 7
## [467,] 5
## [468,] 4
## [469,] 6
## [470,] 5
## [471,] 5
## [472,] 8
## [473,] 6
## [474,] 5
## [475,] 6
## [476,] 7
## [477,] 4
## [478,] 5
## [479,] 7
## [480,] 6
## [481,] 7
## [482,] 5
## [483,] 6
## [484,] 6
## [485,] 5
## [486,] 6
## [487,] 5
## [488,] 4
## [489,] 6
## [490,] 6
## [491,] 8
## [492,] 4
## [493,] 2
## [494,] 3
## [495,] 5
## [496,] 7
## [497,] 6
## [498,] 5
## [499,] 3
## [500,] 3
## [501,] 4
## [502,] 5
## [503,] 7
## [504,] 5
## [505,] 4
## [506,] 6
## [507,] 8
## [508,] 6
## [509,] 6
## [510,] 7
## [511,] 8
## [512,] 6
## [513,] 6
## [514,] 6
## [515,] 6
## [516,] 6
## [517,] 6
## [518,] 7
## [519,] 4
## [520,] 4
## [521,] 2
## [522,] 1
## [523,] 3
## [524,] 5
## [525,] 4
## [526,] 8
## [527,] 8
## [528,] 5
## [529,] 4
## [530,] 7
## [531,] 7
## [532,] 9
## [533,] 8
## [534,] 5
## [535,] 6
## [536,] 5
## [537,] 6
## [538,] 7
## [539,] 6
## [540,] 7
## [541,] 3
## [542,] 5
## [543,] 6
## [544,] 7
## [545,] 7
## [546,] 8
## [547,] 5
## [548,] 6
## [549,] 6
## [550,] 4
## [551,] 8
## [552,] 7
## [553,] 3
## [554,] 6
## [555,] 6
## [556,] 4
## [557,] 3
## [558,] 3
## [559,] 3
## [560,] 4
## [561,] 4
## [562,] 6
## [563,] 6
## [564,] 5
## [565,] 6
## [566,] 4
## [567,] 5
## [568,] 7
## [569,] 8
## [570,] 7
## [571,] 7
## [572,] 2
## [573,] 4
## [574,] 4
## [575,] 5
## [576,] 8
## [577,] 6
## [578,] 7
## [579,] 5
## [580,] 6
## [581,] 5
## [582,] 5
## [583,] 6
## [584,] 6
## [585,] 6
## [586,] 5
## [587,] 2
## [588,] 5
## [589,] 5
## [590,] 2
## [591,] 5
## [592,] 5
## [593,] 7
## [594,] 6
## [595,] 4
## [596,] 6
## [597,] 6
## [598,] 4
## [599,] 5
## [600,] 8
## [601,] 7
## [602,] 4
## [603,] 5
## [604,] 7
## [605,] 6
## [606,] 8
## [607,] 8
## [608,] 2
## [609,] 8
## [610,] 6
## [611,] 2
## [612,] 9
## [613,] 5
## [614,] 5
## [615,] 5
## [616,] 3
## [617,] 4
## [618,] 8
## [619,] 8
## [620,] 7
## [621,] 5
## [622,] 7
## [623,] 6
## [624,] 7
## [625,] 3
## [626,] 4
## [627,] 6
## [628,] 9
## [629,] 4
## [630,] 4
## [631,] 3
## [632,] 7
## [633,] 8
## [634,] 5
## [635,] 5
## [636,] 5
## [637,] 6
## [638,] 0
## [639,] 0
## [640,] 3
## [641,] 5
## [642,] 5
## [643,] 6
## [644,] 7
## [645,] 5
## [646,] 5
## [647,] 3
## [648,] 4
## [649,] 7
## [650,] 5
## [651,] 8
## [652,] 6
## [653,] 7
## [654,] 2
## [655,] 4
## [656,] 7
## [657,] 9
## [658,] 8
## [659,] 6
## [660,] 8
## [661,] 7
## [662,] 5
## [663,] 4
## [664,] 5
## [665,] 8
## [666,] 7
## [667,] 9
## [668,] 8
## [669,] 7
## [670,] 5
## [671,] 7
## [672,] 7
## [673,] 4
## [674,] 4
## [675,] 6
## [676,] 8
## [677,] 7
## [678,] 5
## [679,] 7
## [680,] 3
## [681,] 5
## [682,] 6
## [683,] 5
## [684,] 7
## [685,] 4
## [686,] 7
## [687,] 6
## [688,] 5
## [689,] 8
## [690,] 5
## [691,] 2
## [692,] 2
## [693,] 1
## [694,] 5
## [695,] 6
## [696,] 8
## [697,] 6
## [698,] 2
## [699,] 3
## [700,] 4
## [701,] 6
## [702,] 5
## [703,] 6
## [704,] 3
## [705,] 6
## [706,] 3
## [707,] 5
## [708,] 6
## [709,] 8
## [710,] 8
## [711,] 4
## [712,] 7
## [713,] 7
## [714,] 5
## [715,] 6
## [716,] 7
## [717,] 7
## [718,] 4
## [719,] 7
## [720,] 3
## [721,] 4
## [722,] 2
## [723,] 5
## [724,] 5
## [725,] 6
## [726,] 5
## [727,] 4
## [728,] 4
## [729,] 4
## [730,] 4
## [731,] 3
## [732,] 5
## [733,] 6
## [734,] 5
## [735,] 6
## [736,] 4
## [737,] 6
## [738,] 6
## [739,] 6
## [740,] 8
## [741,] 7
## [742,] 6
## [743,] 9
## [744,] 9
## [745,] 9
## [746,] 8
## [747,] 5
## [748,] 3
## [749,] 8
## [750,] 7
## [751,] 4
## [752,] 6
## [753,] 5
## [754,] 7
## [755,] 5
## [756,] 5
## [757,] 7
## [758,] 5
## [759,] 7
## [760,] 5
## [761,] 6
## [762,] 9
## [763,] 5
## [764,] 8
## [765,] 7
## [766,] 7
## [767,] 2
## [768,] 8
## [769,] 7
## [770,] 4
## [771,] 5
## [772,] 6
## [773,] 8
## [774,] 7
## [775,] 8
## [776,] 7
## [777,] 7
## [778,] 5
## [779,] 5
## [780,] 4
## [781,] 3
## [782,] 7
## [783,] 5
## [784,] 3
## [785,] 4
## [786,] 3
## [787,] 3
## [788,] 5
## [789,] 5
## [790,] 8
## [791,] 9
## [792,] 8
## [793,] 9
## [794,] 8
## [795,] 5
## [796,] 8
## [797,] 8
## [798,] 6
## [799,] 7
## [800,] 6
## [801,] 8
## [802,] 9
## [803,] 9
## [804,] 8
## [805,] 6
## [806,] 8
## [807,] 6
## [808,] 7
## [809,] 8
## [810,] 8
## [811,] 4
## [812,] 6
## [813,] 7
## [814,] 9
## [815,] 7
## [816,] 4
## [817,] 9
## [818,] 6
## [819,] 5
## [820,] 6
## [821,] 3
## [822,] 5
## [823,] 6
## [824,] 7
## [825,] 5
## [826,] 6
## [827,] 3
## [828,] 3
## [829,] 5
## [830,] 5
## [831,] 5
## [832,] 5
## [833,] 8
## [834,] 4
## [835,] 6
## [836,] 4
## [837,] 7
## [838,] 1
## [839,] 4
## [840,] 2
## [841,] 8
## [842,] 8
## [843,] 8
## [844,] 8
## [845,] 6
## [846,] 3
## [847,] 4
## [848,] 3
## [849,] 3
## [850,] 6
## [851,] 8
## [852,] 8
## [853,] 5
## [854,] 7
## [855,] 5
## [856,] 5
## [857,] 4
## [858,] 8
## [859,] 8
## [860,] 8
## [861,] 7
## [862,] 7
## [863,] 5
## [864,] 8
## [865,] 5
## [866,] 7
## [867,] 8
## [868,] 9
## [869,] 6
## [870,] 8
## [871,] 2
## [872,] 7
## [873,] 9
## [874,] 4
## [875,] 3
## [876,] 0
## [877,] 5
## [878,] 4
## [879,] 2
## [880,] 3
## [881,] 0
## [882,] 9
## [883,] 9
## [884,] 7
## [885,] 5
## [886,] 5
## [887,] 5
## [888,] 3
## [889,] 5
## [890,] 6
## [891,] 6
## [892,] 3
## [893,] 7
## [894,] 6
## [895,] 4
## [896,] 7
## [897,] 4
## [898,] 8
## [899,] 6
## [900,] 5
## [901,] 8
## [902,] 6
## [903,] 5
## [904,] 6
## [905,] 5
## [906,] 7
## [907,] 7
## [908,] 9
## [909,] 8
## [910,] 8
## [911,] 7
## [912,] 8
## [913,] 4
## [914,] 9
## [915,] 4
## [916,] 8
## [917,] 5
## [918,] 5
## [919,] 6
## [920,] 4
## [921,] 6
## [922,] 5
## [923,] 6
## [924,] 5
## [925,] 8
## [926,] 6
## [927,] 3
## [928,] 4
## [929,] 5
## [930,] 5
## [931,] 5
## [932,] 7
## [933,] 6
## [934,] 6
## [935,] 2
## [936,] 6
## [937,] 2
## [938,] 7
## [939,] 9
## [940,] 8
## [941,] 2
## [942,] 4
## [943,] 7
## [944,] 8
## [945,] 8
## [946,] 7
## [947,] 4
## [948,] 3
## [949,] 7
## [950,] 4
## [951,] 6
## [952,] 5
## [953,] 3
## [954,] 5
## [955,] 4
## [956,] 2
## [957,] 6
## [958,] 6
## [959,] 8
## [960,] 9
## [961,] 7
## [962,] 7
## [963,] 6
## [964,] 6
## [965,] 3
## [966,] 3
## [967,] 4
## [968,] 4
## [969,] 4
## [970,] 5
## [971,] 6
## [972,] 6
## [973,] 3
## [974,] 6
## [975,] 4
## [976,] 6
## [977,] 6
## [978,] 9
## [979,] 9
## [980,] 9
## [981,] 8
## [982,] 8
## [983,] 4
## [984,] 4
## [985,] 6
## [986,] 9
## [987,] 3
## [988,] 5
## [989,] 2
## [990,] 5
## [991,] 5
## [992,] 6
## [993,] 7
## [994,] 7
## [995,] 7
## [996,] 7
## [997,] 5
## [998,] 1
## [999,] 4
## [1000,] 3
## [1001,] 4
## [1002,] 5
## [1003,] 7
## [1004,] 9
## [1005,] 6
## [1006,] 9
## [1007,] 6
## [1008,] 8
## [1009,] 8
## [1010,] 5
## [1011,] 4
## [1012,] 6
## [1013,] 6
## [1014,] 7
## [1015,] 6
## [1016,] 4
## [1017,] 6
## [1018,] 4
## [1019,] 6
## [1020,] 3
## [1021,] 3
## [1022,] 7
## [1023,] 3
## [1024,] 7
## [1025,] 2
## [1026,] 7
## [1027,] 7
## [1028,] 8
## [1029,] 7
## [1030,] 7
## [1031,] 8
## [1032,] 7
## [1033,] 7
## [1034,] 4
## [1035,] 5
## [1036,] 5
## [1037,] 3
## [1038,] 6
## [1039,] 2
## [1040,] 7
## [1041,] 3
## [1042,] 7
## [1043,] 5
## [1044,] 3
## [1045,] 7
## [1046,] 6
## [1047,] 4
## [1048,] 6
## [1049,] 8
## [1050,] 9
## [1051,] 5
## [1052,] 6
## [1053,] 5
## [1054,] 8
## [1055,] 6
## [1056,] 7
## [1057,] 1
## [1058,] 5
## [1059,] 3
## [1060,] 7
## [1061,] 5
## [1062,] 8
## [1063,] 4
## [1064,] 9
## [1065,] 9
## [1066,] 4
## [1067,] 6
## [1068,] 7
## [1069,] 7
## [1070,] 8
## [1071,] 8
## [1072,] 4
## [1073,] 7
## [1074,] 5
## [1075,] 4
## [1076,] 7
## [1077,] 8
## [1078,] 3
## [1079,] 6
## [1080,] 9
## [1081,] 8
## [1082,] 6
## [1083,] 5
## [1084,] 8
## [1085,] 8
## [1086,] 8
## [1087,] 7
## [1088,] 8
## [1089,] 7
## [1090,] 5
## [1091,] 8
## [1092,] 6
## [1093,] 3
## [1094,] 6
## [1095,] 7
## [1096,] 6
## [1097,] 7
## [1098,] 8
## [1099,] 3
## [1100,] 5
## [1101,] 3
## [1102,] 8
## [1103,] 1
## [1104,] 5
## [1105,] 8
## [1106,] 7
## [1107,] 8
## [1108,] 3
## [1109,] 3
## [1110,] 7
## [1111,] 8
## [1112,] 9
## [1113,] 7
## [1114,] 7
## [1115,] 9
## [1116,] 8
## [1117,] 8
## [1118,] 6
## [1119,] 8
## [1120,] 7
## [1121,] 5
## [1122,] 7
## [1123,] 3
## [1124,] 8
## [1125,] 8
## [1126,] 5
## [1127,] 5
## [1128,] 6
## [1129,] 4
## [1130,] 3
## [1131,] 5
## [1132,] 6
## [1133,] 8
## [1134,] 5
## [1135,] 7
## [1136,] 5
## [1137,] 5
## [1138,] 6
## [1139,] 2
## [1140,] 2
## [1141,] 7
## [1142,] 4
## [1143,] 8
## [1144,] 8
## [1145,] 3
## [1146,] 9
## [1147,] 9
## [1148,] 7
## [1149,] 7
## [1150,] 6
## [1151,] 6
## [1152,] 5
## [1153,] 6
## [1154,] 9
## [1155,] 9
## [1156,] 9
## [1157,] 7
## [1158,] 2
## [1159,] 8
## [1160,] 6
## [1161,] 4
## [1162,] 7
## [1163,] 2
## [1164,] 6
## [1165,] 6
## [1166,] 3
## [1167,] 5
## [1168,] 3
## [1169,] 6
## [1170,] 8
## [1171,] 6
## [1172,] 7
## [1173,] 3
## [1174,] 2
## [1175,] 4
## [1176,] 9
## [1177,] 8
## [1178,] 9
## [1179,] 7
## [1180,] 8
## [1181,] 4
## [1182,] 6
## [1183,] 6
## [1184,] 9
## [1185,] 8
## [1186,] 8
## [1187,] 4
## [1188,] 4
## [1189,] 6
## [1190,] 6
## [1191,] 3
## [1192,] 7
## [1193,] 3
## [1194,] 5
## [1195,] 7
## [1196,] 5
## [1197,] 7
## [1198,] 8
## [1199,] 9
## [1200,] 5
## [1201,] 5
## [1202,] 9
## [1203,] 2
## [1204,] 2
## [1205,] 4
## [1206,] 3
## [1207,] 4
## [1208,] 5
## [1209,] 7
## [1210,] 8
## [1211,] 6
## [1212,] 7
## [1213,] 8
## [1214,] 8
## [1215,] 8
## [1216,] 7
## [1217,] 5
## [1218,] 7
## [1219,] 8
## [1220,] 8
## [1221,] 6
## [1222,] 8
## [1223,] 5
## [1224,] 7
## [1225,] 6
## [1226,] 6
## [1227,] 6
## [1228,] 5
## [1229,] 4
## [1230,] 6
## [1231,] 5
## [1232,] 7
## [1233,] 5
## [1234,] 1
## [1235,] 3
## [1236,] 3
## [1237,] 3
## [1238,] 2
## [1239,] 8
## [1240,] 6
## [1241,] 9
## [1242,] 9
## [1243,] 9
## [1244,] 8
## [1245,] 6
## [1246,] 9
## [1247,] 8
## [1248,] 7
## [1249,] 5
## [1250,] 5
## [1251,] 7
## [1252,] 5
## [1253,] 9
## [1254,] 7
## [1255,] 9
## [1256,] 7
## [1257,] 6
## [1258,] 5
## [1259,] 7
## [1260,] 7
## [1261,] 9
## [1262,] 6
## [1263,] 2
## [1264,] 3
## [1265,] 2
## [1266,] 3
## [1267,] 3
## [1268,] 7
## [1269,] 6
## [1270,] 6
## [1271,] 5
## [1272,] 6
## [1273,] 6
## [1274,] 4
## [1275,] 7
## [1276,] 7
## [1277,] 5
## [1278,] 4
## [1279,] 6
## [1280,] 7
## [1281,] 4
## [1282,] 8
## [1283,] 3
## [1284,] 5
## [1285,] 2
## [1286,] 4
## [1287,] 6
## [1288,] 3
## [1289,] 5
## [1290,] 8
## [1291,] 6
## [1292,] 9
## [1293,] 5
## [1294,] 7
## [1295,] 6
## [1296,] 8
## [1297,] 4
## [1298,] 8
## [1299,] 8
## [1300,] 6
## [1301,] 6
## [1302,] 5
## [1303,] 7
## [1304,] 5
## [1305,] 7
## [1306,] 9
## [1307,] 7
## [1308,] 8
## [1309,] 7
## [1310,] 8
## [1311,] 6
## [1312,] 9
## [1313,] 8
## [1314,] 5
## [1315,] 4
## [1316,] 4
## [1317,] 7
## [1318,] 4
## [1319,] 8
## [1320,] 7
## [1321,] 6
## [1322,] 7
## [1323,] 5
## [1324,] 6
## [1325,] 6
## [1326,] 7
## [1327,] 4
## [1328,] 8
## [1329,] 8
## [1330,] 6
## [1331,] 8
## [1332,] 6
## [1333,] 4
## [1334,] 5
## [1335,] 4
## [1336,] 5
## [1337,] 6
## [1338,] 3
## [1339,] 7
## [1340,] 8
## [1341,] 8
## [1342,] 9
## [1343,] 8
## [1344,] 7
## [1345,] 6
## [1346,] 6
## [1347,] 7
## [1348,] 4
## [1349,] 6
## [1350,] 4
## [1351,] 8
## [1352,] 4
## [1353,] 7
## [1354,] 6
## [1355,] 7
## [1356,] 7
## [1357,] 6
## [1358,] 3
## [1359,] 5
## [1360,] 2
## [1361,] 5
## [1362,] 3
## [1363,] 2
## [1364,] 6
## [1365,] 2
## [1366,] 6
## [1367,] 6
## [1368,] 3
## [1369,] 6
## [1370,] 5
## [1371,] 5
## [1372,] 1
## [1373,] 7
## [1374,] 4
## [1375,] 8
## [1376,] 6
## [1377,] 7
## [1378,] 2
## [1379,] 3
## [1380,] 3
## [1381,] 5
## [1382,] 6
## [1383,] 4
## [1384,] 5
## [1385,] 8
## [1386,] 7
## [1387,] 8
## [1388,] 3
## [1389,] 7
## [1390,] 6
## [1391,] 9
## [1392,] 4
## [1393,] 7
## [1394,] 4
## [1395,] 4
## [1396,] 8
## [1397,] 7
## [1398,] 6
## [1399,] 6
## [1400,] 5
## [1401,] 7
## [1402,] 3
## [1403,] 7
## [1404,] 6
## [1405,] 7
## [1406,] 5
## [1407,] 6
## [1408,] 6
## [1409,] 3
## [1410,] 9
## [1411,] 9
## [1412,] 7
## [1413,] 4
## [1414,] 6
## [1415,] 5
## [1416,] 6
## [1417,] 7
## [1418,] 8
## [1419,] 6
## [1420,] 8
## [1421,] 8
## [1422,] 8
## [1423,] 4
## [1424,] 6
## [1425,] 5
## [1426,] 6
## [1427,] 4
## [1428,] 5
## [1429,] 5
## [1430,] 6
## [1431,] 4
## [1432,] 3
## [1433,] 4
## [1434,] 4
## [1435,] 5
## [1436,] 3
## [1437,] 2
## [1438,] 7
## [1439,] 3
## [1440,] 6
## [1441,] 6
## [1442,] 4
## [1443,] 6
## [1444,] 7
## [1445,] 4
## [1446,] 6
## [1447,] 6
## [1448,] 8
## [1449,] 7
## [1450,] 2
## [1451,] 7
## [1452,] 7
## [1453,] 7
## [1454,] 3
## [1455,] 8
## [1456,] 7
## [1457,] 8
## [1458,] 5
## [1459,] 8
## [1460,] 3
## [1461,] 4
## [1462,] 6
## [1463,] 2
## [1464,] 5
## [1465,] 3
## [1466,] 4
## [1467,] 7
## [1468,] 6
## [1469,] 4
## [1470,] 6
## [1471,] 4
## [1472,] 9
## [1473,] 5
## [1474,] 7
## [1475,] 3
## [1476,] 4
## [1477,] 7
## [1478,] 5
## [1479,] 9
## [1480,] 5
## [1481,] 7
## [1482,] 8
## [1483,] 5
## [1484,] 8
## [1485,] 8
## [1486,] 9
## [1487,] 3
## [1488,] 3
## [1489,] 7
## [1490,] 6
## [1491,] 5
## [1492,] 5
## [1493,] 5
## [1494,] 6
## [1495,] 6
## [1496,] 5
## [1497,] 7
## [1498,] 6
## [1499,] 5
## [1500,] 6
## [1501,] 6
## [1502,] 7
## [1503,] 7
## [1504,] 5
## [1505,] 7
## [1506,] 8
## [1507,] 7
## [1508,] 7
## [1509,] 6
## [1510,] 4
## [1511,] 4
## [1512,] 5
## [1513,] 2
## [1514,] 6
## [1515,] 8
## [1516,] 7
## [1517,] 2
## [1518,] 5
## [1519,] 5
## [1520,] 4
## [1521,] 3
## [1522,] 1
## [1523,] 7
## [1524,] 6
## [1525,] 7
## [1526,] 6
## [1527,] 5
## [1528,] 5
## [1529,] 9
## [1530,] 8
## [1531,] 4
## [1532,] 5
## [1533,] 5
## [1534,] 6
## [1535,] 7
## [1536,] 4
## [1537,] 8
## [1538,] 5
## [1539,] 5
## [1540,] 5
## [1541,] 5
## [1542,] 4
## [1543,] 4
## [1544,] 3
## [1545,] 2
## [1546,] 5
## [1547,] 3
## [1548,] 6
## [1549,] 6
## [1550,] 7
## [1551,] 8
## [1552,] 3
## [1553,] 4
## [1554,] 5
## [1555,] 6
## [1556,] 5
## [1557,] 6
## [1558,] 5
## [1559,] 5
## [1560,] 6
## [1561,] 8
## [1562,] 5
## [1563,] 6
## [1564,] 4
## [1565,] 5
## [1566,] 6
## [1567,] 6
## [1568,] 7
## [1569,] 8
## [1570,] 5
## [1571,] 2
## [1572,] 5
## [1573,] 9
## [1574,] 9
## [1575,] 8
## [1576,] 8
## [1577,] 7
## [1578,] 5
## [1579,] 2
## [1580,] 2
## [1581,] 5
## [1582,] 4
## [1583,] 5
## [1584,] 4
## [1585,] 6
## [1586,] 7
## [1587,] 5
## [1588,] 7
## [1589,] 5
## [1590,] 6
## [1591,] 5
## [1592,] 2
## [1593,] 7
## [1594,] 5
## [1595,] 5
## [1596,] 5
## [1597,] 7
## [1598,] 7
## [1599,] 6
## [1600,] 6
## [1601,] 7
## [1602,] 6
## [1603,] 6
## [1604,] 8
## [1605,] 7
## [1606,] 7
## [1607,] 6
## [1608,] 6
## [1609,] 8
## [1610,] 6
## [1611,] 6
## [1612,] 2
## [1613,] 5
## [1614,] 5
## [1615,] 5
## [1616,] 4
## [1617,] 4
## [1618,] 7
## [1619,] 4
## [1620,] 6
## [1621,] 7
## [1622,] 7
## [1623,] 6
## [1624,] 9
## [1625,] 6
## [1626,] 2
## [1627,] 2
## [1628,] 5
## [1629,] 4
## [1630,] 8
## [1631,] 6
## [1632,] 4
## [1633,] 5
## [1634,] 8
## [1635,] 3
## [1636,] 5
## [1637,] 7
## [1638,] 3
## [1639,] 4
## [1640,] 4
## [1641,] 4
## [1642,] 7
## [1643,] 5
## [1644,] 4
## [1645,] 6
## [1646,] 7
## [1647,] 2
## [1648,] 5
## [1649,] 7
## [1650,] 8
## [1651,] 7
## [1652,] 4
## [1653,] 5
## [1654,] 8
## [1655,] 8
## [1656,] 5
## [1657,] 6
## [1658,] 6
## [1659,] 6
## [1660,] 4
## [1661,] 5
## [1662,] 1
## [1663,] 4
## [1664,] 3
## [1665,] 2
## [1666,] 6
## [1667,] 3
## [1668,] 2
## [1669,] 8
## [1670,] 8
## [1671,] 4
## [1672,] 8
## [1673,] 5
## [1674,] 6
## [1675,] 5
## [1676,] 5
## [1677,] 7
## [1678,] 6
## [1679,] 2
## [1680,] 1
## [1681,] 6
## [1682,] 5
## [1683,] 6
## [1684,] 6
## [1685,] 7
## [1686,] 7
## [1687,] 4
## [1688,] 9
## [1689,] 4
## [1690,] 5
## [1691,] 5
## [1692,] 5
## [1693,] 2
## [1694,] 5
## [1695,] 8
## [1696,] 6
## [1697,] 8
## [1698,] 4
## [1699,] 8
## [1700,] 8
## [1701,] 5
## [1702,] 5
## [1703,] 9
## [1704,] 7
## [1705,] 7
## [1706,] 2
## [1707,] 8
## [1708,] 5
## [1709,] 7
## [1710,] 6
## [1711,] 8
## [1712,] 5
## [1713,] 5
## [1714,] 4
## [1715,] 8
## [1716,] 5
## [1717,] 3
## [1718,] 4
## [1719,] 6
## [1720,] 8
## [1721,] 7
## [1722,] 8
## [1723,] 7
## [1724,] 5
## [1725,] 5
## [1726,] 5
## [1727,] 7
## [1728,] 8
## [1729,] 7
## [1730,] 5
## [1731,] 5
## [1732,] 1
## [1733,] 5
## [1734,] 5
## [1735,] 5
## [1736,] 8
## [1737,] 6
## [1738,] 5
## [1739,] 5
## [1740,] 4
## [1741,] 9
## [1742,] 6
## [1743,] 7
## [1744,] 5
## [1745,] 7
## [1746,] 7
## [1747,] 8
## [1748,] 7
## [1749,] 9
## [1750,] 8
## [1751,] 5
## [1752,] 5
## [1753,] 7
## [1754,] 5
## [1755,] 5
## [1756,] 7
## [1757,] 9
## [1758,] 2
## [1759,] 2
## [1760,] 3
## [1761,] 2
## [1762,] 7
## [1763,] 6
## [1764,] 6
## [1765,] 7
## [1766,] 6
## [1767,] 7
## [1768,] 2
## [1769,] 4
## [1770,] 6
## [1771,] 7
## [1772,] 4
## [1773,] 4
## [1774,] 2
## [1775,] 5
## [1776,] 6
## [1777,] 6
## [1778,] 7
## [1779,] 7
## [1780,] 5
## [1781,] 8
## [1782,] 8
## [1783,] 8
## [1784,] 9
## [1785,] 9
## [1786,] 7
## [1787,] 6
## [1788,] 5
## [1789,] 5
## [1790,] 6
## [1791,] 5
## [1792,] 6
## [1793,] 5
## [1794,] 5
## [1795,] 8
## [1796,] 6
## [1797,] 9
## [1798,] 6
## [1799,] 4
## [1800,] 5
## [1801,] 3
## [1802,] 6
## [1803,] 8
## [1804,] 7
## [1805,] 8
## [1806,] 7
## [1807,] 5
## [1808,] 9
## [1809,] 6
## [1810,] 7
## [1811,] 8
## [1812,] 7
## [1813,] 3
## [1814,] 5
## [1815,] 6
## [1816,] 5
## [1817,] 2
## [1818,] 8
## [1819,] 2
## [1820,] 3
## [1821,] 6
## [1822,] 3
## [1823,] 2
## [1824,] 7
## [1825,] 4
## [1826,] 8
## [1827,] 6
## [1828,] 6
## [1829,] 7
## [1830,] 7
## [1831,] 8
## [1832,] 7
## [1833,] 7
## [1834,] 2
## [1835,] 8
## [1836,] 3
## [1837,] 2
## [1838,] 5
## [1839,] 8
## [1840,] 5
## [1841,] 2
## [1842,] 5
## [1843,] 7
## [1844,] 2
## [1845,] 7
## [1846,] 4
## [1847,] 8
## [1848,] 7
## [1849,] 8
## [1850,] 5
## [1851,] 7
## [1852,] 7
## [1853,] 6
## [1854,] 2
## [1855,] 7
## [1856,] 7
## [1857,] 4
## [1858,] 3
## [1859,] 6
## [1860,] 7
## [1861,] 7
## [1862,] 5
## [1863,] 2
## [1864,] 4
## [1865,] 7
## [1866,] 4
## [1867,] 4
## [1868,] 8
## [1869,] 8
## [1870,] 8
## [1871,] 7
## [1872,] 7
## [1873,] 5
## [1874,] 9
## [1875,] 8
## [1876,] 7
## [1877,] 8
## [1878,] 3
## [1879,] 5
## [1880,] 4
## [1881,] 8
## [1882,] 9
## [1883,] 9
## [1884,] 7
## [1885,] 3
## [1886,] 6
## [1887,] 4
## [1888,] 7
## [1889,] 8
## [1890,] 5
## [1891,] 5
## [1892,] 7
## [1893,] 9
## [1894,] 5
## [1895,] 7
## [1896,] 6
## [1897,] 2
## [1898,] 7
## [1899,] 6
## [1900,] 6
## [1901,] 5
## [1902,] 5
## [1903,] 7
## [1904,] 8
## [1905,] 5
## [1906,] 7
## [1907,] 4
## [1908,] 6
## [1909,] 7
## [1910,] 9
## [1911,] 8
## [1912,] 7
## [1913,] 8
## [1914,] 5
## [1915,] 6
## [1916,] 7
## [1917,] 7
## [1918,] 5
## [1919,] 5
## [1920,] 8
## [1921,] 4
## [1922,] 6
## [1923,] 5
## [1924,] 4
## [1925,] 9
## [1926,] 7
## [1927,] 4
## [1928,] 4
## [1929,] 9
## [1930,] 5
## [1931,] 5
## [1932,] 6
## [1933,] 8
## [1934,] 5
## [1935,] 6
## [1936,] 7
## [1937,] 2
## [1938,] 2
## [1939,] 3
## [1940,] 4
## [1941,] 3
## [1942,] 3
## [1943,] 4
## [1944,] 3
## [1945,] 8
## [1946,] 5
## [1947,] 4
## [1948,] 3
## [1949,] 5
## [1950,] 5
## [1951,] 8
## [1952,] 8
## [1953,] 5
## [1954,] 4
## [1955,] 3
## [1956,] 7
## [1957,] 8
## [1958,] 7
## [1959,] 5
## [1960,] 9
## [1961,] 8
## [1962,] 1
## [1963,] 8
## [1964,] 4
## [1965,] 7
## [1966,] 6
## [1967,] 6
## [1968,] 6
## [1969,] 6
## [1970,] 7
## [1971,] 5
## [1972,] 6
## [1973,] 3
## [1974,] 5
## [1975,] 7
## [1976,] 7
## [1977,] 5
## [1978,] 5
## [1979,] 5
## [1980,] 5
## [1981,] 6
## [1982,] 3
## [1983,] 5
## [1984,] 5
## [1985,] 7
## [1986,] 8
## [1987,] 7
## [1988,] 8
## [1989,] 7
## [1990,] 9
## [1991,] 5
## [1992,] 5
## [1993,] 6
## [1994,] 5
## [1995,] 1
## [1996,] 3
## [1997,] 2
## [1998,] 2
## [1999,] 4
## [2000,] 9
## [2001,] 8
## [2002,] 4
## [2003,] 6
## [2004,] 6
## [2005,] 5
## [2006,] 6
## [2007,] 8
## [2008,] 9
## [2009,] 7
## [2010,] 6
## [2011,] 3
## [2012,] 6
## [2013,] 7
## [2014,] 8
## [2015,] 6
## [2016,] 9
## [2017,] 6
## [2018,] 5
## [2019,] 6
## [2020,] 5
## [2021,] 8
## [2022,] 2
## [2023,] 8
## [2024,] 9
## [2025,] 5
## [2026,] 7
## [2027,] 6
## [2028,] 2
## [2029,] 7
## [2030,] 6
## [2031,] 5
## [2032,] 9
## [2033,] 6
## [2034,] 5
## [2035,] 5
## [2036,] 4
## [2037,] 5
## [2038,] 5
## [2039,] 6
## [2040,] 5
## [2041,] 5
## [2042,] 3
## [2043,] 6
## [2044,] 7
## [2045,] 4
## [2046,] 7
## [2047,] 7
## [2048,] 5
## [2049,] 2
## [2050,] 2
## [2051,] 5
## [2052,] 3
## [2053,] 5
## [2054,] 2
## [2055,] 3
## [2056,] 8
## [2057,] 6
## [2058,] 5
## [2059,] 5
## [2060,] 4
## [2061,] 8
## [2062,] 3
## [2063,] 2
## [2064,] 8
## [2065,] 8
## [2066,] 8
## [2067,] 6
## [2068,] 6
## [2069,] 5
## [2070,] 9
## [2071,] 8
## [2072,] 8
## [2073,] 5
## [2074,] 4
## [2075,] 4
## [2076,] 2
## [2077,] 4
## [2078,] 3
## [2079,] 3
## [2080,] 4
## [2081,] 5
## [2082,] 3
## [2083,] 5
## [2084,] 6
## [2085,] 7
## [2086,] 5
## [2087,] 6
## [2088,] 6
## [2089,] 7
## [2090,] 5
## [2091,] 6
## [2092,] 8
## [2093,] 6
## [2094,] 4
## [2095,] 4
## [2096,] 5
## [2097,] 5
## [2098,] 9
## [2099,] 9
## [2100,] 7
## [2101,] 4
## [2102,] 7
## [2103,] 5
## [2104,] 3
## [2105,] 5
## [2106,] 7
## [2107,] 7
## [2108,] 8
## [2109,] 8
## [2110,] 7
## [2111,] 7
## [2112,] 7
## [2113,] 7
## [2114,] 8
## [2115,] 7
## [2116,] 6
## [2117,] 8
## [2118,] 8
## [2119,] 6
## [2120,] 4
## [2121,] 5
## [2122,] 4
## [2123,] 2
## [2124,] 1
## [2125,] 5
## [2126,] 4
## [2127,] 6
## [2128,] 7
## [2129,] 6
## [2130,] 3
## [2131,] 4
## [2132,] 3
## [2133,] 5
## [2134,] 6
## [2135,] 5
## [2136,] 5
## [2137,] 7
## [2138,] 7
## [2139,] 7
## [2140,] 6
## [2141,] 7
## [2142,] 6
## [2143,] 7
## [2144,] 4
## [2145,] 7
## [2146,] 7
## [2147,] 3
## [2148,] 5
## [2149,] 6
## [2150,] 5
## [2151,] 8
## [2152,] 6
## [2153,] 4
## [2154,] 7
## [2155,] 5
## [2156,] 8
## [2157,] 6
## [2158,] 5
## [2159,] 6
## [2160,] 8
## [2161,] 4
## [2162,] 6
## [2163,] 6
## [2164,] 9
## [2165,] 5
## [2166,] 4
## [2167,] 5
## [2168,] 4
## [2169,] 7
## [2170,] 7
## [2171,] 6
## [2172,] 2
## [2173,] 8
## [2174,] 6
## [2175,] 7
## [2176,] 4
## [2177,] 5
## [2178,] 6
## [2179,] 7
## [2180,] 4
## [2181,] 7
## [2182,] 5
## [2183,] 3
## [2184,] 5
## [2185,] 3
## [2186,] 3
## [2187,] 5
## [2188,] 8
## [2189,] 9
## [2190,] 9
## [2191,] 7
## [2192,] 6
## [2193,] 7
## [2194,] 5
## [2195,] 8
## [2196,] 6
## [2197,] 3
## [2198,] 9
## [2199,] 9
## [2200,] 8
## [2201,] 7
## [2202,] 5
## [2203,] 6
## [2204,] 6
## [2205,] 4
## [2206,] 5
## [2207,] 5
## [2208,] 6
## [2209,] 5
## [2210,] 6
## [2211,] 4
## [2212,] 8
## [2213,] 4
## [2214,] 3
## [2215,] 4
## [2216,] 2
## [2217,] 5
## [2218,] 7
## [2219,] 4
## [2220,] 6
## [2221,] 6
## [2222,] 4
## [2223,] 5
## [2224,] 3
## [2225,] 4
## [2226,] 5
## [2227,] 5
## [2228,] 6
## [2229,] 4
## [2230,] 4
## [2231,] 5
## [2232,] 7
## [2233,] 7
## [2234,] 6
## [2235,] 5
## [2236,] 5
## [2237,] 4
## [2238,] 6
## [2239,] 5
## [2240,] 8
## [2241,] 9
## [2242,] 7
## [2243,] 6
## [2244,] 7
## [2245,] 4
## [2246,] 6
## [2247,] 6
## [2248,] 8
## [2249,] 5
## [2250,] 9
## [2251,] 7
## [2252,] 5
## [2253,] 8
## [2254,] 4
## [2255,] 4
## [2256,] 3
## [2257,] 5
## [2258,] 8
## [2259,] 5
## [2260,] 8
## [2261,] 5
## [2262,] 6
## [2263,] 7
## [2264,] 6
## [2265,] 7
## [2266,] 8
## [2267,] 5
## [2268,] 4
## [2269,] 5
## [2270,] 4
## [2271,] 8
## [2272,] 6
## [2273,] 7
## [2274,] 4
## [2275,] 7
## [2276,] 9
## [2277,] 6
## [2278,] 6
## [2279,] 4
## [2280,] 4
## [2281,] 5
## [2282,] 7
## [2283,] 9
## [2284,] 4
## [2285,] 5
## [2286,] 6
## [2287,] 9
## [2288,] 7
## [2289,] 8
## [2290,] 5
## [2291,] 7
## [2292,] 8
## [2293,] 6
## [2294,] 8
## [2295,] 8
## [2296,] 7
## [2297,] 5
## [2298,] 7
## [2299,] 6
## [2300,] 8
## [2301,] 5
## [2302,] 6
## [2303,] 7
## [2304,] 5
## [2305,] 5
## [2306,] 5
## [2307,] 4
## [2308,] 3
## [2309,] 6
## [2310,] 6
## [2311,] 5
## [2312,] 7
## [2313,] 5
## [2314,] 5
## [2315,] 2
## [2316,] 4
## [2317,] 3
## [2318,] 5
## [2319,] 4
## [2320,] 6
## [2321,] 5
## [2322,] 5
## [2323,] 5
## [2324,] 6
## [2325,] 7
## [2326,] 7
## [2327,] 6
## [2328,] 6
## [2329,] 9
## [2330,] 5
## [2331,] 7
## [2332,] 4
## [2333,] 6
## [2334,] 6
## [2335,] 8
## [2336,] 6
## [2337,] 6
## [2338,] 6
## [2339,] 6
## [2340,] 4
## [2341,] 8
## [2342,] 9
## [2343,] 9
## [2344,] 9
## [2345,] 6
## [2346,] 3
## [2347,] 8
## [2348,] 8
## [2349,] 5
## [2350,] 6
## [2351,] 4
## [2352,] 8
## [2353,] 5
## [2354,] 4
## [2355,] 4
## [2356,] 7
## [2357,] 7
## [2358,] 5
## [2359,] 5
## [2360,] 3
## [2361,] 6
## [2362,] 8
## [2363,] 8
## [2364,] 7
## [2365,] 4
## [2366,] 1
## [2367,] 3
## [2368,] 5
## [2369,] 5
## [2370,] 6
## [2371,] 8
## [2372,] 9
## [2373,] 8
## [2374,] 1
## [2375,] 6
## [2376,] 6
## [2377,] 7
## [2378,] 6
## [2379,] 5
## [2380,] 4
## [2381,] 6
## [2382,] 8
## [2383,] 9
## [2384,] 8
## [2385,] 5
## [2386,] 6
## [2387,] 3
## [2388,] 6
## [2389,] 1
## [2390,] 4
## [2391,] 8
## [2392,] 7
## [2393,] 4
## [2394,] 3
## [2395,] 9
## [2396,] 7
## [2397,] 5
## [2398,] 7
## [2399,] 4
## [2400,] 5
## [2401,] 8
## [2402,] 8
## [2403,] 8
## [2404,] 7
## [2405,] 8
## [2406,] 5
## [2407,] 8
## [2408,] 6
## [2409,] 5
## [2410,] 3
## [2411,] 5
## [2412,] 7
## [2413,] 7
## [2414,] 3
## [2415,] 6
## [2416,] 6
## [2417,] 7
## [2418,] 6
## [2419,] 8
## [2420,] 7
## [2421,] 4
## [2422,] 4
## [2423,] 6
## [2424,] 4
## [2425,] 6
## [2426,] 3
## [2427,] 5
## [2428,] 6
## [2429,] 7
## [2430,] 5
## [2431,] 7
## [2432,] 3
## [2433,] 7
## [2434,] 4
## [2435,] 6
## [2436,] 8
## [2437,] 5
## [2438,] 5
## [2439,] 5
## [2440,] 5
## [2441,] 3
## [2442,] 4
## [2443,] 4
## [2444,] 3
## [2445,] 4
## [2446,] 3
## [2447,] 4
## [2448,] 6
## [2449,] 5
## [2450,] 4
## [2451,] 4
## [2452,] 7
## [2453,] 3
## [2454,] 4
## [2455,] 9
## [2456,] 6
## [2457,] 9
## [2458,] 4
## [2459,] 7
## [2460,] 5
## [2461,] 5
## [2462,] 7
## [2463,] 7
## [2464,] 3
## [2465,] 3
## [2466,] 4
## [2467,] 4
## [2468,] 8
## [2469,] 8
## [2470,] 3
## [2471,] 6
## [2472,] 5
## [2473,] 6
## [2474,] 5
## [2475,] 7
## [2476,] 8
## [2477,] 8
## [2478,] 7
## [2479,] 7
## [2480,] 5
## [2481,] 6
## [2482,] 8
## [2483,] 7
## [2484,] 6
## [2485,] 4
## [2486,] 5
## [2487,] 5
## [2488,] 5
## [2489,] 6
## [2490,] 7
## [2491,] 3
## [2492,] 7
## [2493,] 4
## [2494,] 6
## [2495,] 7
## [2496,] 7
## [2497,] 3
## [2498,] 6
## [2499,] 7
## [2500,] 6
## [2501,] 7
## [2502,] 6
## [2503,] 6
## [2504,] 7
## [2505,] 8
## [2506,] 6
## [2507,] 7
## [2508,] 9
## [2509,] 5
## [2510,] 7
## [2511,] 8
## [2512,] 3
## [2513,] 6
## [2514,] 9
## [2515,] 8
## [2516,] 7
## [2517,] 6
## [2518,] 8
## [2519,] 8
## [2520,] 7
## [2521,] 6
## [2522,] 6
## [2523,] 6
## [2524,] 8
## [2525,] 6
## [2526,] 7
## [2527,] 4
## [2528,] 7
## [2529,] 6
## [2530,] 8
## [2531,] 7
## [2532,] 7
## [2533,] 9
## [2534,] 7
## [2535,] 3
## [2536,] 6
## [2537,] 6
## [2538,] 4
## [2539,] 8
## [2540,] 8
## [2541,] 6
## [2542,] 5
## [2543,] 2
## [2544,] 6
## [2545,] 3
## [2546,] 3
## [2547,] 5
## [2548,] 8
## [2549,] 5
## [2550,] 4
## [2551,] 5
## [2552,] 6
## [2553,] 6
## [2554,] 5
## [2555,] 4
## [2556,] 8
## [2557,] 6
## [2558,] 8
## [2559,] 8
## [2560,] 3
## [2561,] 8
## [2562,] 8
## [2563,] 4
## [2564,] 8
## [2565,] 7
## [2566,] 7
## [2567,] 8
## [2568,] 8
## [2569,] 6
## [2570,] 6
## [2571,] 8
## [2572,] 5
## [2573,] 7
## [2574,] 6
## [2575,] 3
## [2576,] 7
## [2577,] 5
## [2578,] 7
## [2579,] 9
## [2580,] 3
## [2581,] 6
## [2582,] 4
## [2583,] 9
## [2584,] 4
## [2585,] 6
## [2586,] 7
## [2587,] 8
## [2588,] 7
## [2589,] 8
## [2590,] 6
## [2591,] 8
## [2592,] 4
## [2593,] 6
## [2594,] 3
## [2595,] 5
## [2596,] 6
## [2597,] 5
## [2598,] 6
## [2599,] 5
## [2600,] 2
## [2601,] 4
## [2602,] 5
## [2603,] 4
## [2604,] 6
## [2605,] 5
## [2606,] 5
## [2607,] 5
## [2608,] 5
## [2609,] 9
## [2610,] 4
## [2611,] 7
## [2612,] 5
## [2613,] 1
## [2614,] 7
## [2615,] 6
## [2616,] 5
## [2617,] 3
## [2618,] 5
## [2619,] 6
## [2620,] 7
## [2621,] 3
## [2622,] 5
## [2623,] 5
## [2624,] 4
## [2625,] 5
## [2626,] 6
## [2627,] 5
## [2628,] 6
## [2629,] 8
## [2630,] 8
## [2631,] 5
## [2632,] 4
## [2633,] 4
## [2634,] 4
## [2635,] 4
## [2636,] 6
## [2637,] 7
## [2638,] 4
## [2639,] 7
## [2640,] 6
## [2641,] 7
## [2642,] 8
## [2643,] 8
## [2644,] 8
## [2645,] 7
## [2646,] 7
## [2647,] 4
## [2648,] 4
## [2649,] 9
## [2650,] 7
## [2651,] 3
## [2652,] 2
## [2653,] 2
## [2654,] 3
## [2655,] 4
## [2656,] 6
## [2657,] 5
## [2658,] 2
## [2659,] 2
## [2660,] 7
## [2661,] 5
## [2662,] 7
## [2663,] 5
## [2664,] 6
## [2665,] 6
## [2666,] 6
## [2667,] 8
## [2668,] 7
## [2669,] 8
## [2670,] 8
## [2671,] 4
## [2672,] 8
## [2673,] 4
## [2674,] 5
## [2675,] 7
## [2676,] 7
## [2677,] 6
## [2678,] 5
## [2679,] 7
## [2680,] 6
## [2681,] 6
## [2682,] 8
## [2683,] 4
## [2684,] 7
## [2685,] 7
## [2686,] 6
## [2687,] 6
## [2688,] 5
## [2689,] 9
## [2690,] 7
## [2691,] 5
## [2692,] 7
## [2693,] 7
## [2694,] 8
## [2695,] 5
## [2696,] 6
## [2697,] 3
## [2698,] 1
## [2699,] 5
## [2700,] 4
## [2701,] 5
## [2702,] 2
## [2703,] 4
## [2704,] 5
## [2705,] 2
## [2706,] 1
## [2707,] 3
## [2708,] 6
## [2709,] 4
## [2710,] 8
## [2711,] 5
## [2712,] 4
## [2713,] 5
## [2714,] 5
## [2715,] 8
## [2716,] 6
## [2717,] 6
## [2718,] 5
## [2719,] 4
## [2720,] 5
## [2721,] 7
## [2722,] 3
## [2723,] 3
## [2724,] 4
## [2725,] 8
## [2726,] 8
## [2727,] 6
## [2728,] 7
## [2729,] 5
## [2730,] 5
## [2731,] 6
## [2732,] 6
## [2733,] 7
## [2734,] 7
## [2735,] 5
## [2736,] 7
## [2737,] 5
## [2738,] 5
## [2739,] 7
## [2740,] 6
## [2741,] 4
## [2742,] 4
## [2743,] 3
## [2744,] 4
## [2745,] 2
## [2746,] 8
## [2747,] 6
## [2748,] 6
## [2749,] 8
## [2750,] 6
## [2751,] 8
## [2752,] 5
## [2753,] 3
## [2754,] 5
## [2755,] 5
## [2756,] 7
## [2757,] 9
## [2758,] 5
## [2759,] 7
## [2760,] 5
## [2761,] 8
## [2762,] 9
## [2763,] 8
## [2764,] 7
## [2765,] 8
## [2766,] 6
## [2767,] 5
## [2768,] 6
## [2769,] 4
## [2770,] 4
## [2771,] 6
## [2772,] 5
## [2773,] 7
## [2774,] 5
## [2775,] 3
## [2776,] 9
## [2777,] 8
## [2778,] 5
## [2779,] 8
## [2780,] 6
## [2781,] 8
## [2782,] 5
## [2783,] 6
## [2784,] 6
## [2785,] 3
## [2786,] 6
## [2787,] 5
## [2788,] 6
## [2789,] 4
## [2790,] 8
## [2791,] 4
## [2792,] 3
## [2793,] 5
## [2794,] 1
## [2795,] 4
## [2796,] 3
## [2797,] 7
## [2798,] 6
## [2799,] 6
## [2800,] 7
## [2801,] 5
## [2802,] 6
## [2803,] 5
## [2804,] 3
## [2805,] 5
## [2806,] 7
## [2807,] 8
## [2808,] 7
## [2809,] 6
## [2810,] 5
## [2811,] 2
## [2812,] 5
## [2813,] 6
## [2814,] 6
## [2815,] 9
## [2816,] 7
## [2817,] 8
## [2818,] 7
## [2819,] 5
## [2820,] 7
## [2821,] 8
## [2822,] 8
## [2823,] 7
## [2824,] 5
## [2825,] 6
## [2826,] 6
## [2827,] 1
## [2828,] 4
## [2829,] 3
## [2830,] 4
## [2831,] 7
## [2832,] 7
## [2833,] 7
## [2834,] 6
## [2835,] 7
## [2836,] 8
## [2837,] 7
## [2838,] 6
## [2839,] 6
## [2840,] 5
## [2841,] 5
## [2842,] 5
## [2843,] 6
## [2844,] 6
## [2845,] 6
## [2846,] 5
## [2847,] 4
## [2848,] 4
## [2849,] 5
## [2850,] 6
## [2851,] 5
## [2852,] 6
## [2853,] 5
## [2854,] 7
## [2855,] 4
## [2856,] 5
## [2857,] 4
## [2858,] 2
## [2859,] 3
## [2860,] 6
## [2861,] 5
## [2862,] 3
## [2863,] 5
## [2864,] 3
## [2865,] 6
## [2866,] 3
## [2867,] 4
## [2868,] 8
## [2869,] 7
## [2870,] 3
## [2871,] 5
## [2872,] 5
## [2873,] 7
## [2874,] 6
## [2875,] 2
## [2876,] 2
## [2877,] 3
## [2878,] 3
## [2879,] 8
## [2880,] 8
## [2881,] 6
## [2882,] 4
## [2883,] 7
## [2884,] 8
## [2885,] 7
## [2886,] 7
## [2887,] 6
## [2888,] 4
## [2889,] 6
## [2890,] 6
## [2891,] 8
## [2892,] 3
## [2893,] 3
## [2894,] 8
## [2895,] 6
## [2896,] 6
## [2897,] 5
## [2898,] 3
## [2899,] 5
## [2900,] 4
## [2901,] 8
## [2902,] 7
## [2903,] 5
## [2904,] 6
## [2905,] 9
## [2906,] 7
## [2907,] 7
## [2908,] 5
## [2909,] 5
## [2910,] 3
## [2911,] 3
## [2912,] 0
## [2913,] 5
## [2914,] 1
## [2915,] 3
## [2916,] 8
## [2917,] 8
## [2918,] 5
## [2919,] 6
## [2920,] 5
## [2921,] 7
## [2922,] 7
## [2923,] 7
## [2924,] 8
## [2925,] 6
## [2926,] 5
## [2927,] 3
## [2928,] 5
## [2929,] 8
## [2930,] 4
## [2931,] 6
## [2932,] 7
## [2933,] 4
## [2934,] 5
## [2935,] 6
## [2936,] 3
## [2937,] 5
## [2938,] 3
## [2939,] 6
## [2940,] 7
## [2941,] 8
## [2942,] 6
## [2943,] 4
## [2944,] 5
## [2945,] 6
## [2946,] 5
## [2947,] 3
## [2948,] 8
## [2949,] 6
## [2950,] 6
## [2951,] 6
## [2952,] 5
## [2953,] 7
## [2954,] 7
## [2955,] 8
## [2956,] 6
## [2957,] 7
## [2958,] 3
## [2959,] 5
## [2960,] 3
## [2961,] 6
## [2962,] 5
## [2963,] 3
## [2964,] 7
## [2965,] 8
## [2966,] 7
## [2967,] 9
## [2968,] 9
## [2969,] 8
## [2970,] 8
## [2971,] 5
## [2972,] 4
## [2973,] 2
## [2974,] 7
## [2975,] 7
## [2976,] 8
## [2977,] 4
## [2978,] 7
## [2979,] 7
## [2980,] 1
## [2981,] 7
## [2982,] 2
## [2983,] 4
## [2984,] 2
## [2985,] 8
## [2986,] 5
## [2987,] 6
## [2988,] 9
## [2989,] 4
## [2990,] 7
## [2991,] 5
## [2992,] 5
## [2993,] 2
## [2994,] 6
## [2995,] 4
## [2996,] 1
## [2997,] 3
## [2998,] 6
## [2999,] 6
## [3000,] 4
## [3001,] 4
## [3002,] 5
## [3003,] 8
## [3004,] 7
## [3005,] 8
## [3006,] 3
## [3007,] 4
## [3008,] 4
## [3009,] 5
## [3010,] 6
## [3011,] 3
## [3012,] 3
## [3013,] 9
## [3014,] 6
## [3015,] 7
## [3016,] 7
## [3017,] 5
## [3018,] 7
## [3019,] 5
## [3020,] 6
## [3021,] 3
## [3022,] 7
## [3023,] 5
## [3024,] 3
## [3025,] 4
## [3026,] 4
## [3027,] 3
## [3028,] 7
## [3029,] 7
## [3030,] 7
## [3031,] 8
## [3032,] 7
## [3033,] 8
## [3034,] 5
## [3035,] 2
## [3036,] 3
## [3037,] 6
## [3038,] 8
## [3039,] 6
## [3040,] 4
## [3041,] 7
## [3042,] 3
## [3043,] 5
## [3044,] 6
## [3045,] 5
## [3046,] 7
## [3047,] 6
## [3048,] 8
## [3049,] 6
## [3050,] 7
## [3051,] 6
## [3052,] 5
## [3053,] 5
## [3054,] 4
## [3055,] 4
## [3056,] 7
## [3057,] 6
## [3058,] 7
## [3059,] 1
## [3060,] 1
## [3061,] 6
## [3062,] 6
## [3063,] 6
## [3064,] 3
## [3065,] 6
## [3066,] 6
## [3067,] 7
## [3068,] 7
## [3069,] 3
## [3070,] 6
## [3071,] 7
## [3072,] 8
## [3073,] 7
## [3074,] 7
## [3075,] 9
## [3076,] 7
## [3077,] 6
## [3078,] 6
## [3079,] 9
## [3080,] 7
## [3081,] 7
## [3082,] 4
## [3083,] 4
## [3084,] 5
## [3085,] 2
## [3086,] 7
## [3087,] 4
## [3088,] 6
## [3089,] 8
## [3090,] 9
## [3091,] 6
## [3092,] 5
## [3093,] 7
## [3094,] 4
## [3095,] 5
## [3096,] 6
## [3097,] 4
## [3098,] 8
## [3099,] 7
## [3100,] 3
## [3101,] 3
## [3102,] 1
## [3103,] 9
## [3104,] 7
## [3105,] 7
## [3106,] 4
## [3107,] 7
## [3108,] 4
## [3109,] 7
## [3110,] 5
## [3111,] 5
## [3112,] 7
## [3113,] 6
## [3114,] 9
## [3115,] 5
## [3116,] 5
## [3117,] 5
## [3118,] 6
## [3119,] 7
## [3120,] 8
## [3121,] 4
## [3122,] 9
## [3123,] 6
## [3124,] 3
## [3125,] 2
## [3126,] 6
## [3127,] 6
## [3128,] 7
## [3129,] 4
## [3130,] 4
## [3131,] 5
## [3132,] 3
## [3133,] 6
## [3134,] 8
## [3135,] 4
## [3136,] 7
## [3137,] 5
## [3138,] 5
## [3139,] 5
## [3140,] 6
## [3141,] 5
## [3142,] 7
## [3143,] 4
## [3144,] 1
## [3145,] 4
## [3146,] 4
## [3147,] 5
## [3148,] 5
## [3149,] 8
## [3150,] 8
## [3151,] 5
## [3152,] 5
## [3153,] 5
## [3154,] 3
## [3155,] 6
## [3156,] 5
## [3157,] 4
## [3158,] 3
## [3159,] 4
## [3160,] 3
## [3161,] 8
## [3162,] 4
## [3163,] 9
## [3164,] 3
## [3165,] 6
## [3166,] 9
## [3167,] 6
## [3168,] 5
## [3169,] 6
## [3170,] 8
## [3171,] 7
## [3172,] 5
## [3173,] 8
## [3174,] 8
## [3175,] 4
## [3176,] 7
## [3177,] 6
## [3178,] 2
## [3179,] 2
## [3180,] 6
## [3181,] 6
## [3182,] 8
## [3183,] 7
## [3184,] 6
## [3185,] 7
## [3186,] 5
## [3187,] 8
## [3188,] 7
## [3189,] 7
## [3190,] 7
## [3191,] 4
## [3192,] 4
## [3193,] 5
## [3194,] 3
## [3195,] 4
## [3196,] 7
## [3197,] 6
## [3198,] 5
## [3199,] 6
## [3200,] 8
## [3201,] 4
## [3202,] 7
## [3203,] 9
## [3204,] 7
## [3205,] 6
## [3206,] 7
## [3207,] 5
## [3208,] 7
## [3209,] 5
## [3210,] 2
## [3211,] 6
## [3212,] 3
## [3213,] 6
## [3214,] 5
## [3215,] 7
## [3216,] 6
## [3217,] 6
## [3218,] 5
## [3219,] 7
## [3220,] 4
## [3221,] 3
## [3222,] 6
## [3223,] 5
## [3224,] 2
## [3225,] 5
## [3226,] 7
## [3227,] 6
## [3228,] 2
## [3229,] 7
## [3230,] 4
## [3231,] 3
## [3232,] 2
## [3233,] 3
## [3234,] 8
## [3235,] 6
## [3236,] 9
## [3237,] 3
## [3238,] 8
## [3239,] 6
## [3240,] 4
## [3241,] 4
## [3242,] 6
## [3243,] 6
## [3244,] 7
## [3245,] 8
## [3246,] 2
## [3247,] 9
## [3248,] 7
## [3249,] 6
## [3250,] 6
## [3251,] 4
## [3252,] 8
## [3253,] 8
## [3254,] 7
## [3255,] 4
## [3256,] 7
## [3257,] 6
## [3258,] 7
## [3259,] 5
## [3260,] 8
## [3261,] 6
## [3262,] 6
## [3263,] 6
## [3264,] 9
## [3265,] 6
## [3266,] 6
## [3267,] 7
## [3268,] 7
## [3269,] 5
## [3270,] 4
## [3271,] 4
## [3272,] 5
## [3273,] 5
## [3274,] 6
## [3275,] 9
## [3276,] 8
## [3277,] 7
## [3278,] 3
## [3279,] 1
## [3280,] 9
## [3281,] 8
## [3282,] 8
## [3283,] 8
## [3284,] 7
## [3285,] 7
## [3286,] 6
## [3287,] 2
## [3288,] 5
## [3289,] 5
## [3290,] 5
## [3291,] 7
## [3292,] 7
## [3293,] 8
## [3294,] 4
## [3295,] 5
## [3296,] 5
## [3297,] 5
## [3298,] 6
## [3299,] 6
## [3300,] 6
## [3301,] 7
## [3302,] 8
## [3303,] 6
## [3304,] 8
## [3305,] 5
## [3306,] 4
## [3307,] 4
## [3308,] 7
## [3309,] 7
## [3310,] 4
## [3311,] 4
## [3312,] 5
## [3313,] 4
## [3314,] 4
## [3315,] 4
## [3316,] 4
## [3317,] 7
## [3318,] 6
## [3319,] 5
## [3320,] 7
## [3321,] 6
## [3322,] 3
## [3323,] 8
## [3324,] 4
## [3325,] 7
## [3326,] 5
## [3327,] 5
## [3328,] 5
## [3329,] 5
## [3330,] 5
## [3331,] 5
## [3332,] 7
## [3333,] 6
## [3334,] 4
## [3335,] 8
## [3336,] 7
## [3337,] 7
## [3338,] 6
## [3339,] 7
## [3340,] 6
## [3341,] 8
## [3342,] 5
## [3343,] 6
## [3344,] 7
## [3345,] 5
## [3346,] 4
## [3347,] 1
## [3348,] 8
## [3349,] 8
## [3350,] 5
## [3351,] 4
## [3352,] 5
## [3353,] 5
## [3354,] 3
## [3355,] 1
## [3356,] 6
## [3357,] 5
## [3358,] 6
## [3359,] 6
## [3360,] 5
## [3361,] 9
## [3362,] 4
## [3363,] 5
## [3364,] 7
## [3365,] 4
## [3366,] 6
## [3367,] 6
## [3368,] 7
## [3369,] 6
## [3370,] 6
## [3371,] 7
## [3372,] 9
## [3373,] 9
## [3374,] 8
## [3375,] 8
## [3376,] 9
## [3377,] 4
## [3378,] 6
## [3379,] 5
## [3380,] 6
## [3381,] 8
## [3382,] 5
## [3383,] 4
## [3384,] 7
## [3385,] 5
## [3386,] 3
## [3387,] 5
## [3388,] 4
## [3389,] 5
## [3390,] 3
## [3391,] 8
## [3392,] 7
## [3393,] 3
## [3394,] 8
## [3395,] 8
## [3396,] 8
## [3397,] 7
## [3398,] 9
## [3399,] 8
## [3400,] 8
## [3401,] 5
## [3402,] 6
## [3403,] 8
## [3404,] 8
## [3405,] 7
## [3406,] 9
## [3407,] 5
## [3408,] 5
## [3409,] 7
## [3410,] 5
## [3411,] 7
## [3412,] 5
## [3413,] 6
## [3414,] 4
## [3415,] 2
## [3416,] 3
## [3417,] 3
## [3418,] 9
## [3419,] 4
## [3420,] 7
## [3421,] 8
## [3422,] 6
## [3423,] 7
## [3424,] 7
## [3425,] 7
## [3426,] 8
## [3427,] 7
## [3428,] 9
## [3429,] 7
## [3430,] 9
## [3431,] 9
## [3432,] 7
## [3433,] 9
## [3434,] 7
## [3435,] 8
## [3436,] 8
## [3437,] 4
## [3438,] 5
## [3439,] 5
## [3440,] 4
## [3441,] 6
## [3442,] 5
## [3443,] 4
## [3444,] 7
## [3445,] 6
## [3446,] 7
## [3447,] 7
## [3448,] 1
## [3449,] 5
## [3450,] 5
## [3451,] 6
## [3452,] 7
## [3453,] 7
## [3454,] 7
## [3455,] 8
## [3456,] 8
## [3457,] 8
## [3458,] 7
## [3459,] 7
## [3460,] 5
## [3461,] 9
## [3462,] 7
## [3463,] 9
## [3464,] 1
## [3465,] 2
## [3466,] 3
## [3467,] 4
## [3468,] 2
## [3469,] 5
## [3470,] 3
## [3471,] 3
## [3472,] 3
## [3473,] 5
## [3474,] 5
## [3475,] 7
## [3476,] 8
## [3477,] 8
## [3478,] 8
## [3479,] 3
## [3480,] 5
## [3481,] 5
## [3482,] 6
## [3483,] 8
## [3484,] 5
## [3485,] 3
## [3486,] 8
## [3487,] 6
## [3488,] 7
## [3489,] 7
## [3490,] 6
## [3491,] 9
## [3492,] 2
## [3493,] 5
## [3494,] 4
## [3495,] 6
## [3496,] 5
## [3497,] 5
## [3498,] 3
## [3499,] 3
## [3500,] 7
## [3501,] 7
## [3502,] 7
## [3503,] 7
## [3504,] 6
## [3505,] 9
## [3506,] 6
## [3507,] 6
## [3508,] 4
## [3509,] 8
## [3510,] 5
## [3511,] 9
## [3512,] 5
## [3513,] 8
## [3514,] 7
## [3515,] 4
## [3516,] 7
## [3517,] 5
## [3518,] 4
## [3519,] 4
## [3520,] 4
## [3521,] 7
## [3522,] 5
## [3523,] 7
## [3524,] 8
## [3525,] 8
## [3526,] 7
## [3527,] 2
## [3528,] 5
## [3529,] 4
## [3530,] 6
## [3531,] 5
## [3532,] 2
## [3533,] 4
## [3534,] 7
## [3535,] 9
## [3536,] 6
## [3537,] 8
## [3538,] 6
## [3539,] 6
## [3540,] 8
## [3541,] 5
## [3542,] 8
## [3543,] 3
## [3544,] 5
## [3545,] 6
## [3546,] 2
## [3547,] 5
## [3548,] 8
## [3549,] 7
## [3550,] 8
## [3551,] 9
## [3552,] 8
## [3553,] 7
## [3554,] 7
## [3555,] 7
## [3556,] 5
## [3557,] 7
## [3558,] 6
## [3559,] 8
## [3560,] 6
## [3561,] 8
## [3562,] 8
## [3563,] 9
## [3564,] 7
## [3565,] 8
## [3566,] 5
## [3567,] 5
## [3568,] 5
## [3569,] 4
## [3570,] 7
## [3571,] 6
## [3572,] 7
## [3573,] 4
## [3574,] 8
## [3575,] 8
## [3576,] 7
## [3577,] 7
## [3578,] 3
## [3579,] 7
## [3580,] 8
## [3581,] 8
## [3582,] 5
## [3583,] 8
## [3584,] 8
## [3585,] 4
## [3586,] 8
## [3587,] 5
## [3588,] 6
## [3589,] 6
## [3590,] 3
## [3591,] 5
## [3592,] 4
## [3593,] 7
## [3594,] 5
## [3595,] 2
## [3596,] 3
## [3597,] 7
## [3598,] 8
## [3599,] 7
## [3600,] 3
## [3601,] 9
## [3602,] 4
## [3603,] 6
## [3604,] 6
## [3605,] 8
## [3606,] 8
## [3607,] 2
## [3608,] 4
## [3609,] 6
## [3610,] 6
## [3611,] 2
## [3612,] 8
## [3613,] 4
## [3614,] 5
## [3615,] 5
## [3616,] 1
## [3617,] 7
## [3618,] 5
## [3619,] 4
## [3620,] 6
## [3621,] 7
## [3622,] 7
## [3623,] 5
## [3624,] 4
## [3625,] 8
## [3626,] 5
## [3627,] 7
## [3628,] 9
## [3629,] 8
## [3630,] 8
## [3631,] 3
## [3632,] 6
## [3633,] 7
## [3634,] 6
## [3635,] 4
## [3636,] 5
## [3637,] 5
## [3638,] 6
## [3639,] 4
## [3640,] 6
## [3641,] 4
## [3642,] 7
## [3643,] 5
## [3644,] 5
## [3645,] 6
## [3646,] 6
## [3647,] 6
## [3648,] 7
## [3649,] 7
## [3650,] 7
## [3651,] 7
## [3652,] 9
## [3653,] 7
## [3654,] 8
## [3655,] 8
## [3656,] 8
## [3657,] 8
## [3658,] 8
## [3659,] 5
## [3660,] 5
## [3661,] 6
## [3662,] 3
## [3663,] 8
## [3664,] 3
## [3665,] 4
## [3666,] 6
## [3667,] 5
## [3668,] 7
## [3669,] 7
## [3670,] 8
## [3671,] 5
## [3672,] 4
## [3673,] 4
## [3674,] 5
## [3675,] 8
## [3676,] 3
## [3677,] 6
## [3678,] 9
## [3679,] 5
## [3680,] 1
## [3681,] 9
## [3682,] 7
## [3683,] 6
## [3684,] 4
## [3685,] 7
## [3686,] 6
## [3687,] 8
## [3688,] 5
## [3689,] 6
## [3690,] 5
## [3691,] 9
## [3692,] 5
## [3693,] 6
## [3694,] 7
## [3695,] 7
## [3696,] 8
## [3697,] 6
## [3698,] 4
## [3699,] 6
## [3700,] 7
## [3701,] 6
## [3702,] 4
## [3703,] 5
## [3704,] 5
## [3705,] 5
## [3706,] 6
## [3707,] 7
## [3708,] 3
## [3709,] 6
## [3710,] 8
## [3711,] 9
## [3712,] 9
## [3713,] 0
## [3714,] 2
## [3715,] 8
## [3716,] 6
## [3717,] 7
## [3718,] 9
## [3719,] 6
## [3720,] 7
## [3721,] 9
## [3722,] 5
## [3723,] 7
## [3724,] 5
## [3725,] 8
## [3726,] 4
## [3727,] 5
## [3728,] 4
## [3729,] 8
## [3730,] 9
## [3731,] 7
## [3732,] 6
## [3733,] 8
## [3734,] 9
## [3735,] 4
## [3736,] 5
## [3737,] 6
## [3738,] 7
## [3739,] 6
## [3740,] 4
## [3741,] 5
## [3742,] 6
## [3743,] 6
## [3744,] 3
## [3745,] 2
## [3746,] 6
## [3747,] 5
## [3748,] 7
## [3749,] 6
## [3750,] 9
## [3751,] 9
## [3752,] 7
## [3753,] 6
## [3754,] 7
## [3755,] 8
## [3756,] 7
## [3757,] 9
## [3758,] 7
## [3759,] 7
## [3760,] 7
## [3761,] 6
## [3762,] 6
## [3763,] 8
## [3764,] 5
## [3765,] 7
## [3766,] 5
## [3767,] 7
## [3768,] 6
## [3769,] 4
## [3770,] 4
## [3771,] 8
## [3772,] 7
## [3773,] 4
## [3774,] 5
## [3775,] 7
## [3776,] 9
## [3777,] 6
## [3778,] 6
## [3779,] 9
## [3780,] 9
## [3781,] 6
## [3782,] 7
## [3783,] 3
## [3784,] 4
## [3785,] 2
## [3786,] 2
## [3787,] 4
## [3788,] 4
## [3789,] 4
## [3790,] 5
## [3791,] 3
## [3792,] 5
## [3793,] 6
## [3794,] 6
## [3795,] 6
## [3796,] 5
## [3797,] 5
## [3798,] 6
## [3799,] 5
## [3800,] 4
## [3801,] 8
## [3802,] 8
## [3803,] 8
## [3804,] 7
## [3805,] 6
## [3806,] 3
## [3807,] 6
## [3808,] 2
## [3809,] 7
## [3810,] 3
## [3811,] 3
## [3812,] 4
## [3813,] 6
## [3814,] 7
## [3815,] 6
## [3816,] 3
## [3817,] 6
## [3818,] 5
## [3819,] 5
## [3820,] 3
## [3821,] 5
## [3822,] 8
## [3823,] 0
## [3824,] 7
## [3825,] 7
## [3826,] 8
## [3827,] 4
## [3828,] 6
## [3829,] 5
## [3830,] 3
## [3831,] 3
## [3832,] 3
## [3833,] 2
## [3834,] 7
## [3835,] 7
## [3836,] 7
## [3837,] 5
## [3838,] 8
## [3839,] 6
## [3840,] 5
## [3841,] 3
## [3842,] 7
## [3843,] 6
## [3844,] 4
## [3845,] 6
## [3846,] 8
## [3847,] 6
## [3848,] 7
## [3849,] 7
## [3850,] 5
## [3851,] 5
## [3852,] 2
## [3853,] 5
## [3854,] 5
## [3855,] 6
## [3856,] 8
## [3857,] 8
## [3858,] 8
## [3859,] 9
## [3860,] 8
## [3861,] 7
## [3862,] 9
## [3863,] 4
## [3864,] 7
## [3865,] 7
## [3866,] 7
## [3867,] 6
## [3868,] 8
## [3869,] 8
## [3870,] 6
## [3871,] 7
## [3872,] 8
## [3873,] 6
## [3874,] 8
## [3875,] 5
## [3876,] 5
## [3877,] 7
## [3878,] 9
## [3879,] 6
## [3880,] 5
## [3881,] 6
## [3882,] 6
## [3883,] 4
## [3884,] 2
## [3885,] 6
## [3886,] 9
## [3887,] 6
## [3888,] 8
## [3889,] 8
## [3890,] 6
## [3891,] 7
## [3892,] 8
## [3893,] 6
## [3894,] 6
## [3895,] 7
## [3896,] 5
## [3897,] 3
## [3898,] 7
## [3899,] 7
## [3900,] 7
## [3901,] 9
## [3902,] 6
## [3903,] 6
## [3904,] 7
## [3905,] 7
## [3906,] 6
## [3907,] 2
## [3908,] 7
## [3909,] 5
## [3910,] 3
## [3911,] 7
## [3912,] 6
## [3913,] 7
## [3914,] 6
## [3915,] 7
## [3916,] 7
## [3917,] 6
## [3918,] 8
## [3919,] 7
## [3920,] 5
## [3921,] 7
## [3922,] 7
## [3923,] 5
## [3924,] 7
## [3925,] 8
## [3926,] 6
## [3927,] 4
## [3928,] 6
## [3929,] 7
## [3930,] 4
## [3931,] 4
## [3932,] 3
## [3933,] 7
## [3934,] 6
## [3935,] 8
## [3936,] 5
## [3937,] 4
## [3938,] 6
## [3939,] 7
## [3940,] 5
## [3941,] 6
## [3942,] 3
## [3943,] 5
## [3944,] 3
## [3945,] 5
## [3946,] 7
## [3947,] 4
## [3948,] 8
## [3949,] 2
## [3950,] 4
## [3951,] 7
## [3952,] 8
## [3953,] 7
## [3954,] 5
## [3955,] 3
## [3956,] 6
## [3957,] 7
## [3958,] 4
## [3959,] 3
## [3960,] 8
## [3961,] 5
## [3962,] 7
## [3963,] 7
## [3964,] 8
## [3965,] 6
## [3966,] 5
## [3967,] 3
## [3968,] 9
## [3969,] 5
## [3970,] 2
## [3971,] 2
## [3972,] 4
## [3973,] 1
## [3974,] 4
## [3975,] 5
## [3976,] 8
## [3977,] 4
## [3978,] 4
## [3979,] 5
## [3980,] 9
## [3981,] 4
## [3982,] 5
## [3983,] 5
## [3984,] 7
## [3985,] 4
## [3986,] 9
## [3987,] 5
## [3988,] 7
## [3989,] 7
## [3990,] 4
## [3991,] 7
## [3992,] 0
## [3993,] 4
## [3994,] 5
## [3995,] 5
## [3996,] 7
## [3997,] 6
## [3998,] 9
## [3999,] 8
## [4000,] 7
data.frame(obs = ppd) %>%
ggplot(aes(x=obs)) +
geom_histogram() +
labs(x="Observed water (out of 9)")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Oops, we’ve jumped ahead of ourselves! Best practice is to simulate values from your prior first and check to see if those priors are reasonable.
m1p <-
brm(data = list(w = 6),
family = binomial(link = "identity"),
w | trials(9) ~ 0 + Intercept,
prior(beta(1, 1), class = b, lb = 0, ub = 1),
iter = 5000, warmup = 1000, seed = 3, chains=1,
sample_prior = "only")
samples_from_prior = as_draws_df(m1p)
samples_from_prior %>%
ggplot(aes(x=b_Intercept)) +
geom_density(fill = "grey", color = "white") +
labs(x="Proportion water", title="Prior")
We may also want the PRIOR PREDICTIVE DISTRIBUTION which is the expected observiations given our prior.
prior_pd = posterior_predict(m1p)
data.frame(obs = prior_pd) %>%
ggplot(aes(x=obs)) +
geom_histogram() +
labs(x="Observed water (out of 9)")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Simulating from your priors – prior predictive simulation – is an essential part of modeling. This allows you to see what your choices imply about the data. You’ll be able to diagnose bad choices.
At this point in the course, I’m going to start throwing a lot of code at you. Do I expect you to memorize this code? Of course not.
Do you need to understand every single thing that’s happening in the code? Nope.
But, you’ll learn a lot by taking the time to figure out what’s happening in a code chunk. Class time will frequently include exercises where I ask you to adapt code I’ve shared in the slides to a new dataset or to answer a new problem. When doing so, go back through the old code and figure out what’s going on. Run the code one line at a time. Always observe the output and take some time to look at the object that was created or modified. Here are some functions that will be extremely useful:
str() # what kind of object is this? what is its structure?
dim() # what are the dimensions (rows/columns) of this object
head() # give me the first bit of this object
str(prior_pd)
## int [1:4000, 1] 6 0 1 1 9 7 0 5 4 5 ...
## - attr(*, "dimnames")=List of 2
## ..$ : NULL
## ..$ : NULL
dim(prior_pd)
## [1] 4000 1
head(prior_pd)
## [,1]
## [1,] 6
## [2,] 0
## [3,] 1
## [4,] 1
## [5,] 9
## [6,] 7
The globe tossing example is cute and easy to work with, but let’s move towards the kinds of variables we more frequently work with. Let’s create a model for some outcome, \(y\) that is continuous.
\[\begin{align*} y_i &\sim \text{Normal}(\mu, \sigma) \\ \mu &\sim \text{Normal}(0, 10) \\ \sigma &\sim \text{Uniform}(0, 5) \end{align*}\]
set.seed(9)
y = rnorm(n = 31, mean = 4, sd = .5)
m2 = brm(
data = list(y=y),
family = gaussian,
y ~ 1,
prior = c(prior( normal(0,10), class=Intercept),
prior( uniform(0,5), class=sigma)),
iter = 5000, warmup = 1000, seed = 3, chains=1,
file = here("files/models/m21.2")
)
Using the Howell data (don’t load the rethinking package
because it can interfere with brms).
data("Howell1", package = "rethinking")
d <- Howell1
str(d)
## 'data.frame': 544 obs. of 4 variables:
## $ height: num 152 140 137 157 145 ...
## $ weight: num 47.8 36.5 31.9 53 41.3 ...
## $ age : num 63 63 65 41 51 35 32 27 19 54 ...
## $ male : int 1 0 0 1 0 1 0 1 0 1 ...
library(measurements)
d$height <- conv_unit(d$height, from = "cm", to = "feet")
d$weight <- conv_unit(d$weight, from = "kg", to = "lbs")
rethinking::precis(d)
## mean sd 5.5% 94.5% histogram
## height 4.5362072 0.9055921 2.661042 5.4375 ▁▁▁▁▂▂▇▇▁
## weight 78.5079631 32.4502290 20.636856 120.1583 ▁▂▃▂▂▁▁▃▅▇▇▃▂▁
## age 29.3443934 20.7468882 1.000000 66.1350 ▇▅▅▃▅▂▂▁▁
## male 0.4724265 0.4996986 0.000000 1.0000 ▇▁▁▁▁▁▁▁▁▇
d2 <- d[ d$age >= 18, ]
Write a mathematical model for the weights in this data set. (Don’t worry about predicting from other variables yet.)
\[\begin{align*} w &\sim \text{Normal}(\mu, \sigma) \\ \mu &\sim \text{Normal}(130, 20) \\ \sigma &\sim \text{Uniform}(0, 25) \\ \end{align*}\]
\[\begin{align*} w &\sim \text{Normal}(\mu, \sigma) \\ \mu &\sim \text{Normal}(130, 20) \\ \sigma &\sim \text{Uniform}(0, 25) \\ \end{align*}\]
Simulate from your priors (parameters values and prior predictive values).
Sample from your priors:
m3p = brm(
data = d2,
family = gaussian,
weight ~ 1,
prior = c(prior( normal(130,20), class=Intercept),
prior( uniform(0,25), class=sigma, lb=0, ub=25)),
iter = 5000, warmup = 1000, seed = 3, chains=1,
sample_prior = "only")
Sampling parameter estimates for the Intercept.
pairs(m3p)
## Warning: Only one chain in 'x'. This plot is more useful with multiple chains.
This is a different (shorter) way to plot your posterior. Good things: it automatically includes the scatterplot so you can see the implications of how these parameters correlate. Bad things: not customizable and not useable when you have a lot of parameters.
Simulate values of weight.
prior_pd = posterior_predict(m3p)
dim(prior_pd)
## [1] 4000 352
as.data.frame(prior_pd) %>%
pivot_longer(everything()) %>%
ggplot(aes(x=value)) +
geom_histogram() +
labs(x="Expected observed weights (based on prior)")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Another shorter way:
pp_check(m3p)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
m3 = brm(
data = d2,
family = gaussian,
weight ~ 1,
prior = c(prior( normal(130,20), class=Intercept),
prior( uniform(0,25), class=sigma, lb=0, ub=25)),
iter = 5000, warmup = 1000, seed = 3, chains=1,
file = here("files/models/m21.3"))
posterior_summary(m3)
## Estimate Est.Error Q2.5 Q97.5
## b_Intercept 99.221909 0.75720356 97.751074 100.737488
## sigma 14.291987 0.55198322 13.254363 15.428183
## Intercept 99.221909 0.75720356 97.751074 100.737488
## lprior -8.318377 0.05827127 -8.433538 -8.203915
## lp__ -1441.294276 1.02695551 -1444.235605 -1440.292326
pairs(m3)
## Warning: Only one chain in 'x'. This plot is more useful with multiple chains.
posterior_predict(m3) %>%
as.data.frame() %>%
pivot_longer(everything()) %>%
ggplot(aes(x=value)) +
geom_density(fill = "grey", color = "white") +
geom_density( aes(x = weight), data=d2, inherit.aes = F)
pp_check(m3)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
We might assume that height and weight are associated with each other. Indeed, within our sample:
plot(d2$weight ~ d2$height)
Update your mathematical model to incorporate height.
\[\begin{align*} w_i &\sim \text{Normal}(\mu_i, \sigma) \\ \mu_i &= \alpha + \beta h_i \\ \alpha &\sim \text{Normal}(??, ??) \\ \beta &\sim \text{Normal}(0, 25) \\ \sigma &\sim \text{Uniform}(0, 25) \\ \end{align*}\]
\(=\) is deterministic – once we know other variables, \(\mu_i\) is known with certainty
made-up parameters are the targets of learning
Update your mathematical model to incorporate height.
\[\begin{align*} w_i &\sim \text{Normal}(\mu_i, \sigma) \\ \mu_i &= \alpha + \beta (h_i - \bar{h}) \\ \alpha &\sim \text{Normal}(130, 20) \\ \beta &\sim \text{Normal}(0, 25) \\ \sigma &\sim \text{Uniform}(0, 25) \\ \end{align*}\]
\(=\) is deterministic – once we know other variables, \(\mu_i\) is known with certainty
made-up parameters are the targets of learning
\[\begin{align*} w_i &\sim \text{Normal}(\mu_i, \sigma) \\ \mu_i &= \alpha + \beta (h_i - \bar{h}) \\ \alpha &\sim \text{Normal}(130, 20) \\ \beta &\sim \text{Normal}(0, 25) \\ \sigma &\sim \text{Uniform}(0, 25) \\ \end{align*}\]
To update our brms code:
d2$height_c = d2$height - mean(d2$height)
m4p = brm(
data = d2,
family = gaussian,
weight ~ 1 + height_c,
prior = c(prior( normal(130,20), class=Intercept),
prior( normal(0,25), class=b),
prior( uniform(0,25), class=sigma, lb=0, ub=25)),
iter = 5000, warmup = 1000, seed = 3, chains=1,
sample_prior = "only")
set.seed(9)
samples_from_prior = as_draws_df(m4p)
str(samples_from_prior)
## draws_df [4,000 × 9] (S3: draws_df/draws/tbl_df/tbl/data.frame)
## $ b_Intercept: num [1:4000] 167 155.6 149.6 88.7 176.4 ...
## $ b_height_c : num [1:4000] -57.9 -57.2 -18.5 52 -47.3 ...
## $ sigma : num [1:4000] 14.1 11.8 13.3 2.1 21.1 ...
## $ Intercept : num [1:4000] 167 155.6 149.6 88.7 176.4 ...
## $ lprior : num [1:4000] -15.7 -14.7 -12 -15.6 -15.8 ...
## $ lp__ : num [1:4000] -13.8 -12.9 -10.2 -14.9 -14.6 ...
## $ .chain : int [1:4000] 1 1 1 1 1 1 1 1 1 1 ...
## $ .iteration : int [1:4000] 1 2 3 4 5 6 7 8 9 10 ...
## $ .draw : int [1:4000] 1 2 3 4 5 6 7 8 9 10 ...
d2 %>%
ggplot(aes(x=height_c, y=weight)) +
geom_blank() +
geom_abline( aes(intercept=b_Intercept, slope=b_height_c),
data=samples_from_prior[1:50, ],
alpha=.3) +
scale_x_continuous(name = "height(feet)",
breaks=seq(4,6,by=.5)-mean(d2$height),
labels=seq(4,6,by=.5))
Describe in words what’s wrong with our priors.
Slope should not be negative. How can we fix this?
Could use a uniform distribution bounded by 0.
m4p = brm(
data = d2,
family = gaussian,
weight ~ 1 + height_c,
prior = c(prior( normal(130,20), class=Intercept),
prior( uniform(0,25), class=b),
prior( uniform(0,25), class=sigma, lb=0, ub=25)),
iter = 5000, warmup = 1000, seed = 3, chains=1,
sample_prior = "only")
Fit the new weight model to the data.
m4 = brm(
data = d2,
family = gaussian,
weight ~ 1 + height_c,
prior = c(prior( normal(130,20), class=Intercept),
prior( uniform(0,25), class=b),
prior( uniform(0,25), class=sigma, lb=0, ub=25)),
iter = 5000, warmup = 1000, seed = 3, chains=1,
file=here("files/models/m21.4"))
posterior_summary(m4)
## Estimate Est.Error Q2.5 Q97.5
## b_Intercept 99.21508 0.5154139 98.17317 100.20955
## b_height_c 24.74661 0.2609169 24.05150 24.99355
## sigma 10.36090 0.3898531 9.65234 11.12228
## Intercept 99.21508 0.5154139 98.17317 100.20955
## lprior -11.53739 0.0396881 -11.61861 -11.46176
## lp__ -1332.15882 1.3936194 -1335.90639 -1330.52002
Draw lines from the posterior distribution and plot with the data.
set.seed(9)
samples_from_posterior = as_draws_df(m4)
d2%>%
ggplot(aes(x=height_c, y=weight)) +
geom_point(size=.5) +
geom_abline( aes(intercept=b_Intercept, slope=b_height_c),
data=samples_from_posterior[1:50, ],
alpha=.3,
color="#1c5253") +
scale_x_continuous(name = "height(feet)",
breaks=seq(4,6,by=.5)-mean(d2$height),
labels=seq(4,6,by=.5))
A side note: a major concern or critique of Bayesian analysis is that the subjectivity of the priors allow for nefarious behavior. “Putting our thumbs on the scale,” so to speak. But priors are quickly overwhelmed by data. Case in point:
m4e = brm(
data = d2,
family = gaussian,
weight ~ 1 + height_c,
prior = c(prior( normal(130,20), class=Intercept),
prior( normal(-5,5), class=b),
prior( uniform(0,25), class=sigma, lb=0, ub=25)),
iter = 5000, warmup = 1000, seed = 3, chains=1,
file=here("files/models/m21.4e"))
posterior_summary(m4e)
## Estimate Est.Error Q2.5 Q97.5
## b_Intercept 99.197733 0.5035653 98.214092 100.15932
## b_height_c 35.775818 1.8832462 32.177572 39.50070
## sigma 9.529429 0.3687178 8.839379 10.27233
## Intercept 99.197733 0.5035653 98.214092 100.15932
## lprior -44.172476 3.0738968 -50.456409 -38.49994
## lp__ -1334.629214 1.1849503 -1337.645046 -1333.23509
You’ll only really get into trouble with uniform priors that have a boundary, if true population parameter is outside your boundary. A good rule of thumb is to avoid the uniform distribution. We’ll cover other options for priors for \(\sigma\) in future lectures, but as a preview, the exponential distribution works very well for this!